#include "spectrumanalyzer.h"
#include <QDebug>
#include <cmath>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

SpectrumAnalyzer::SpectrumAnalyzer(int numBands, QObject *parent)
    : QObject(parent)
    , m_numBands(numBands)
    , m_fftSize(2048) // Power of 2 for FFT
{
    m_spectrum.resize(m_numBands);
    m_spectrum.fill(0.0f);
}

void SpectrumAnalyzer::processAudioData(const QVector<float> &audioData, int sampleRate)
{
    if (audioData.size() < m_fftSize) {
        return; // Not enough data for FFT
    }

    // Take the latest chunk of audio data
    QVector<float> chunk(m_fftSize);
    int startIdx = audioData.size() - m_fftSize;
    for (int i = 0; i < m_fftSize; i++) {
        chunk[i] = audioData[startIdx + i];
    }

    // Apply Hann window
    QVector<float> windowed = applyHannWindow(chunk);

    // Perform FFT
    QVector<std::complex<float>> fftOutput;
    performFFT(windowed, fftOutput);

    // Calculate spectrum
    calculateSpectrum(fftOutput, sampleRate);

    // Debug: Check spectrum values
    static int debugCount = 0;
    if (debugCount++ % 30 == 0) { // Log every 30 frames (~1 second)
        float maxSpec = 0.0f;
        for (float val : m_spectrum) {
            maxSpec = qMax(maxSpec, val);
        }
        qDebug() << "Spectrum calculated: max value =" << maxSpec;
    }

    emit spectrumReady(m_spectrum);
}

QVector<float> SpectrumAnalyzer::applyHannWindow(const QVector<float> &input)
{
    QVector<float> output(input.size());
    for (int i = 0; i < input.size(); i++) {
        float window = 0.5f * (1.0f - std::cos(2.0f * M_PI * i / (input.size() - 1)));
        output[i] = input[i] * window;
    }
    return output;
}

void SpectrumAnalyzer::performFFT(const QVector<float> &input, 
                                  QVector<std::complex<float>> &output)
{
    int n = input.size();
    output.resize(n);

    // Convert input to complex
    for (int i = 0; i < n; i++) {
        output[i] = std::complex<float>(input[i], 0.0f);
    }

    // Cooley-Tukey FFT algorithm (Radix-2 DIT)
    // Bit-reversal permutation
    int bits = 0;
    int temp = n;
    while (temp > 1) {
        bits++;
        temp >>= 1;
    }

    for (int i = 0; i < n; i++) {
        int j = 0;
        for (int k = 0; k < bits; k++) {
            if (i & (1 << k)) {
                j |= 1 << (bits - 1 - k);
            }
        }
        if (j > i) {
            std::swap(output[i], output[j]);
        }
    }

    // FFT computation
    for (int s = 1; s <= bits; s++) {
        int m = 1 << s;
        int m2 = m >> 1;
        std::complex<float> w(1.0f, 0.0f);
        std::complex<float> wm = std::exp(std::complex<float>(0.0f, -M_PI / m2));

        for (int j = 0; j < m2; j++) {
            for (int k = j; k < n; k += m) {
                std::complex<float> t = w * output[k + m2];
                std::complex<float> u = output[k];
                output[k] = u + t;
                output[k + m2] = u - t;
            }
            w *= wm;
        }
    }
}

void SpectrumAnalyzer::calculateSpectrum(const QVector<std::complex<float>> &fftData, 
                                         int sampleRate)
{
    // We only use the first half of FFT output (positive frequencies)
    int halfSize = fftData.size() / 2;
    
    // Maximum frequency is Nyquist frequency (sampleRate / 2)
    // We want 0 Hz to 20 kHz
    float maxFreq = std::min(20000.0f, sampleRate / 2.0f);
    
    // Divide frequency range into bands
    // Using logarithmic scale for better perceptual distribution
    QVector<float> bandLimits(m_numBands + 1);
    float minFreq = 20.0f; // Start from 20 Hz
    float logMin = std::log10(minFreq);
    float logMax = std::log10(maxFreq);
    
    for (int i = 0; i <= m_numBands; i++) {
        float logFreq = logMin + (logMax - logMin) * i / m_numBands;
        bandLimits[i] = std::pow(10.0f, logFreq);
    }

    // Calculate magnitude for each band
    for (int band = 0; band < m_numBands; band++) {
        float startFreq = bandLimits[band];
        float endFreq = bandLimits[band + 1];
        
        int startBin = static_cast<int>(startFreq * fftData.size() / sampleRate);
        int endBin = static_cast<int>(endFreq * fftData.size() / sampleRate);
        
        startBin = qMax(0, qMin(startBin, halfSize - 1));
        endBin = qMax(startBin + 1, qMin(endBin, halfSize));
        
        float sum = 0.0f;
        int count = 0;
        
        for (int i = startBin; i < endBin; i++) {
            float magnitude = std::abs(fftData[i]);
            sum += magnitude;
            count++;
        }
        
        if (count > 0) {
            float avgMagnitude = sum / count;
            // Convert to dB scale and normalize
            float db = 20.0f * std::log10(avgMagnitude + 1e-10f);
            // Normalize to 0-1 range (assuming -60 dB to 0 dB range)
            m_spectrum[band] = qMax(0.0f, qMin(1.0f, (db + 60.0f) / 60.0f));
        } else {
            m_spectrum[band] = 0.0f;
        }
    }
}